Fluctuations and correlations in sandpile models 
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We perform numerical simulations of the sandpile model for non- vanishing driving fields h and 
dissipation rates e. Unlike simulations performed in the slow driving limit, the unique time scale 
present in our system allows us to measure unambiguously response and correlation functions. We 
discuss the dynamic scaling of the model and show that fluctuation-dissipation relations are not 
obeyed in this system. 
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The sandpile automaton [JD is one of the simplest 
model of avalanche transport, a phenomenon of grow- 
ing experimental and theoretical interest. In the model 
introduced by Bak, Tang and Wiesenfeld (BTW) 0, 
grains of "energy" are injected into the system. Open 
boundary conditions or bulk dissipation insure a bal- 
ance between input and output flow and allow for a non- 
equilibrium stationary state. In the limit of slow external 
driving and small dissipation, which corresponds to an in- 
finite time scale separation between driving and response, 
the model displays an highly fluctuating avalanche be- 
havior, indicative of a critical point. Despite the impres- 
sive theoretical effort devoted to understanding the crit- 
ical behavior of the model , several important issues 
still remain to be addressed. 

Numerical simulations are usually performed under 
slow driving and boundary dissipation, since the limit 
of infinite time scale separation is easily implemented in 
the computer and provides a simple way to access the 
avalanche critical behavior |^ ||]. However, due to the 
presence of two infinitely separated time scales, an un- 
ambiguous definition of dynamic response and correla- 
tion functions is not possible [jl0|]. This hinders a clear 
characterization of the non-equilibrium stationary state 
in terms of static and dynamic response and correlation 
functions. Evaluation of these quantities helps to elu- 
cidate the nature of the critical point and provides a 
test of fluctuation-dissipation relations, at least in some 
weaker sense. Recently, it has been proposed to interpret 
the behavior of sandpile models in analogy with other 
non-equilibrium critical phenomena, such as absorbing 
pha se transitions ]Tl[ | , driven interfaces in random media 
|l2| [lif and branching processes These theoretical 

studies suggest new ways to perform numerical simula- 
tions in which a unique time scale is considered Jll],[l4|,[l6| . 

In this letter, we present numerical simulations of the 
sandpile model for different driving rates h and study how 
the system approaches the critical point when h — > 0. In 
this way, we are able to measure quantities that are not 
accessible in the time scale separation regime. The local 
density of active sites, that can be identified as the order 
parameter of the model |ll| , is homogeneous only in the 



case of bulk dissipation. For boundary dissipation, it dis- 
plays a marked curvature, that was anticipated in Refs. 

and could explain several scaling anomalies found 
in the BTW model. The energy landscape is instead ho- 
mogeneous in both cases and its statistical properties do 
not depend on the dissipation rate e in the limit h — -> 0. 

We measure correlation and response functions in time 
and space domains and observe the scaling of the re- 
lated characteristic lengths and times. We find two 
different characteristic times, implying that fluctuation- 
dissipation relations are not obeyed. We observe, how- 
ever, a well defined scaling behavior and the value of the 
critical exponents are in agreement with recent large scale 
numerical simulations of slowly driven sandpiles 
Finally, the present numerical analysis opens the road 
to future studies to resolve some longstanding problems 
such as the precise identification of universality classes 
for these models || . 

In sandpile models |l[], each site i of a d— dimensional 
lattice bears an integer variable Zi > 0, which we call 
energy. At each time step an energy grain is added on 
a randomly chosen site. When a site reaches or exceeds 
a threshold z c it topples: Zi — > zi — z c , and Zj — ► Zj + 1 
at each of the g nearest neighbors (nn) of i. Each top- 
pling can trigger nn to topple and so on, generating an 
avalanche. The original BTW model is conservative and 
energy is dissipated only at the boundary, i.e. energy 
grains from toppling boundary sites flow out of the sys- 
tem. Infinitely slow driving is implicitly built into the 
model: during the avalanche the energy input stops, until 
the system is again quiescent (no active sites are present), 
so that we can identify two distinct time scales T<j and 
T a , for driving and activity, respectively. A single driv- 
ing time step can in principle be followed by an infinite 
number of avalanche time steps and T a /Td — > 0. For this 
reason, there are two possible definitions for the correla- 
tion function, depending on the choice of the scale used 
to measure time (slow or fast) 10,18]. 



Here we simulate the BTW sandpile model for a non- 
vanishing driving field: each site has a probability h per 
unit time to receive an energy grain, also if active sites 
are present in the system. This defines a unique time 
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step for both driving and activity updating. The param- 
eter h sets the driving rate, and in the limit h — > + we 
recover the slow driving limit; i.e. during the evolution 
of an avalanche the system does not receive energy. We 
consider two possible mechanisms for energy dissipation: 
(i) usual boundary dissipation and (ii) bulk dissipation, 
simulated introducing the probability a that a toppling 
site looses its energy without transferring it to the neigh- 
bors, which corresponds to an effective average dissipa- 
tion e = az c . In case (ii), we impose periodic boundary 
conditions. We use two dimensional lattices with linear 
sizes ranging from L — 64 to L = 300, and parameters in 
the range lO" 6 < h< 10~ 3 and 10~ 3 < e < 10~ 2 . 

The order parameter in sandpile models is the density 
1 [M_of active sites, whose energy is larger than z c 
l| jl2|J14|JT^ |. The dependence of the order parameter 
on the control parameters h and e is readily obtained 
by means of conservation arguments ]T^ |: since energy 
is conserved in the stationary state, the incoming energy 
flux (Ji n ) = hL d must be equal to the dissipated energy 
(Jout) — e (Pa)L d - By equating the two fluxes we obtain 
(Pa) = h/e. In systems with boundary dissipation, the 
effective dissipation rate scales with the system size as 
e ~ L _/x , with p = 2 @, yielding (p a ) ~ hL 2 . It has to 
be noticed that the model is critical just in the double 
limit e — > and h/e — ► 0. The onset of a nonvanish- 
ing driving thus destroys criticality in that it enforces a 
nonzero dissipation. For feCe the cutoff length scaling is 
dominated only by dissipation, while for greater driving 
fields more complicate scaling behaviors occur . 

We study the behavior of the density of active sites 
in the system and measure the stationary average den- 
sity of local energies; i.e. the density pi of sites with 
i energy grains. In Fig. 1 we report the behavior of 
the densities as a function of h/e. For small values of 
h/e we find (p t ) = pf + cji/e + 0((h/e) 2 ), where p\ are 
the values extrapolated from the limit h — > + and are 
given by p° = 0.075(1), p\ = 0.176(1), p\ = 0.307(1) and 
p§ = 0.442(1). These values are in excellent agreement 
with the exact values obtained for the slowly driven sand- 
pile (with boundary dissipation) Q and are independent 
of the dissipation rate. This implies that the energy sub- 
strate over which avalanches propagate is the same in 
the case of bulk and boundary dissipation. For i > 3 we 
obtain p° = and for small h we observe (p a ) ~ (/04), 
while for larger h higher energy levels become populated 
and (p a ) has non negligible contributions coming from 
(pi) with i > 4. Finally we confirm that (p a ) = h/e for 
the whole range of parameters. In the case of boundary 
dissipation we recover (p a ) ~ hL 2 . 

To elucidate the differences between bulk and bound- 
ary dissipation, we measure the local density of active 
sites (p a (r)). In the case of bulk dissipation the density 
profile is flat (p a (f)) = (p a ), while in the case of bound- 
ary dissipation we obtain a surface that can be well ap- 
proximated by a paraboloid (see Fig. 2). This is due to 
the highly inhomogeneous dissipation which imposes a 
zero density of active sites on the lattice boundary and 



corresponds to an elastic interface pinned at the bound- 
aries as discussed in Ref.s jll 14 . This effect can explain 



the anomalies encountered in the numerical evaluation of 
avalanche exponents [|||] and the persistent deviations 
from simple scaling observed in BTW model |l7j] . 

In order to obtain a quantitative description of the 
stationary state, we study the effect on the stationary 
density of a small perturbation in the driving field 

Ap a (r, t) = [ Xh.Ar -r',t - t')Ah(r',t')dr'dt' (1) 



where Xh.e(f, t) is the local response function. In the limit 
h — > + the integrated susceptibility \ = J dtd d rxh.e (r, t) 
scales as the average avalanche size x ~ (s) and the time 
integrated response function scales as [16| 
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,-r/Z 
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where £ is the characteristic length. Since x = dp a /dh 
and p a = h/e, the response function diverges in the limit 
of vanishing driving and dissipation as x — l/ e - By not- 
ing that x ~ £ 2 ; we obtain that £ ~ e _I/ with v = 1/2 
p6| . These results can be obtained in mean-field (MF) 
theory but hold in all dimensions due to conservation 

To measure the response function, we drive the sys- 
tem in the stationary state with a given h and we then 
add n ener gy grains (i.e. Ah = n/L 2 ) on a given 
lattice site |2Cfl. The time integrated response func- 
tion is equivalent to the average difference of activity 
Xh,e( r ) = Ap a {r) = (p a (r))h+Ah - (p a (r))h, where r de- 
notes the distance from the perturbed site. We observe 
that this function decays exponentially as predicted by 
Eq. and measure the correlation length £ (see Fig. 3). 
In the case of bulk dissipation, for small driving fields 
the £ depends only on the dissipation rate and scales as 
£ - e~ v , with v = 0.50 ± 0.01 (see Fig. 4). In the case 
of boundary dissipation, to evaluate Xh,i( r ) we have to 
consider explicitly the spatial inhomogeneity of the sta- 
tionary density: (p a (r)) ^ h/e. We observe also in this 
case that the integrated response function decays expo- 
nentially and defines a correlation length increasing lin- 
early with the lattice size; i.e £ ~ L. This result does not 
agree with the anomalous scaling found in a continuous 
energy sandpile p8| . We perform analogous simulations 
in d = 3 and find that Eq. || is still verified pl[ . 

Furthermore, we study the response function in the 
time domain defined as Xh.e(t) = J drxh,e( r ^) after a 
small variation Ah of the driving field. Also in this case 
we obtain a clear exponential behavior defining the char- 
acteristic time scale r. For small driving field, r scales 
as a function of the dissipation rate as r ~ e~ A , with 
A = 0.75 ± 0.05 (see Fig. 4). We then evaluate the dy- 
namical exponent z = A/v = 1.5 ± 0.1 relating time 
and spatial characteristic length: r ~ £ z . In the limit 
h — ► + , we expect that the critical exponents v and 
z express the divergence of avalanche characteristic size 
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and time, respectively. The numerical results confirm 
this observation |J. For increasing values of h, the driv- 
ing field enters the scaling form and the results will be 
reported elsewhere El]]. 

We now turn to the analysis of the correlation function 
defined as C(r,t) = (p a (r,t)p a (0,0)) — (p a ) 2 ■ In previous 
simulations, performed in the slow driving limit, corre- 
lation functions were usually measured with respect to 
the slow time scale |l| Jl8| , P3| and the fast time scale was 
explored studying the avalanche propagation. The intro- 
duction of non vanishing driving and dissipation allows 
us to bridge the gap between the two regimes. We study 
the correlation function in time and space domains and 
find an exponential decay at long times and distances 
p2| , defining the correlation lengths £ c and t c for space 
and time, respectively. The scaling of these correlation 
lengths is in agreement with the one obtained analyzing 
the response functions (i.e. £ c ~ e _I/ , with v ~ 0.5 and 
t c ~ e~ A , with A ~ 0.75) and confirms the existence of 
a unique critical behavior in time and space (see Fig. 4). 

In order to clarify the interplay between slow and 
fast dynamical modes, we analyze fluctuation-dissipation 
relations. In equilibrium phenomena, the fluctuation- 
dissipation theorem ensures that the response of the sys- 
tem to a small perturbation is related to the correlation 
function. In particular, the response function is given by 



x(t) 



1 dC(t) 
T dt ' 



(3) 



where T is the temperature. Eq. || is strictly verified only 
in equilibrium systems, but it has been recently general- 
ized to some classes of non equilibrium systems, namely 
systems displaying "aging" | j24| . In those examples the 
fluctuation-dissipation relation provides an information 
on an effective non equilibrium temperature that rules 
the dynamical evolution of the system. 

We test Eq. ^ and we find that the usual linear be- 
havior does not hold. On the contrary, we show that 
the parametric plot of versus C(t) defines a power 
law behavior, as shown in the double logarithmic plot 
of Fig. 5. This is striking evidence that the fluctuation- 
dissipation relation does not hold in these systems. Since 
we are in presence of two exponential functions, the linear 
behavior on the logarithmic scale is the signature of two 
different values for characteristic times r c and r for the 
correlation and response function respectively. The slope 
indicates the ratio among the two time scales is given by 
t c /t ~ 0.4 and does not depend on driving and dissi- 
pation rates. This observation reflects the fact that the 
correlation and the response times scale with the same ex- 
ponents with respect to dissipation and define unambigu- 
ously the critical behavior of the model. In particular, it 
implies that the dynamical exponent z ~ 1.5 is unique 
and can be estimated either by measuring avalanche dis- 
tributions or the correlation functions. Previous simula- 
tions |l],[l8| revealed two different dynamical exponents 
in the fast and in the slow time scale. These differences 



are probably due to the ambiguous definition of time in 
the infinite time scale separation limit. 

Finally, we note that it is not possible to define an ef- 
fective temperature for the dynamics of sandpile models. 
It is interesting to compare this observation with a re- 
cent work [^| showing that the stationary state of non- 
equilibrium threshold models, similar to the one stud- 
ied here, can be described by Boltzmann statistics in 
the mean-field limit. The validity of claim of Ref. |2^] 
has been debated in the literature |2qj . We measure 
fluctuation-dissipation relations in a random neighbor 
sandpile model, which is described by mean-field the- 
ory, and find that fluctuation-dissipation relations are 
not satisfied pl[ , in disagreement with the conclusions 
of Ref. |§. 
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FIG. 2. Local density of active sites (p a (x,y)), in the case 
of boundary dissipation; the linear lattice size is L = 100, and 
h= HT 4 . 
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FIG. 1. Mean densities (p;) of sites with i energy grains, vs 
h/e; Inset: mean density (p a ) of active sites, vs h/e. Values of 
h range from 10 -6 to 10~ 3 , and e is between 10 -3 and 10 -2 . 
We observe that (p a ) = h/e, and that the various (pi) depend 
on h and e only through the ratio h/e (see text). 
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FIG. 3. Time integrated response function Xft^o.f to a con- 
stant perturbation as a function of r; the linear lattice size is 
L — 200. The lines are exponential fits. 
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FIG. 4. Characteristic length £ and characteristic time 
lim h ^ + r(h,e) estimated from the spatial response and cor- 
relation functions with bulk dissipation and for various sys- 
tem sizes and dissipation rates. For small dissipation rates 
£ — > co, and larger lattice sizes must be used. The straight 
lines represent the best fits with slope v = 0.5 and A = 0.75 
for £ and r respectively. 
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FIG. 5. Double logarithmic plot of x(t) vs C(t) for various 
values of h and e; the slope of the straight lines is t c /t ~ 0.4. 
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